##################
#   MRP--example using death penalty responses
#################



load("Death penalty megapoll.RData") ##megapoll of values
load("Weights and state level covariates for MRP.RData") ##state level covariates and weights



library(stringr)
##weights.for.post is the matrix of population characteristics
weights.for.post<-as.data.frame(weights.for.post)
weights.for.post$stfip<-str_pad(as.character(weights.for.post$stfip),2,"left","0")

##state covariates is the matrix with state level effects
state.covariates$stfip<-str_pad(state.covariates$stfip,2,"left","0")
load("State crosswalk.RData")
colnames(state.covariates)[1]<-"fipsst"
colnames(weights.for.post)[2]<-"fipsst"


state.covariates<-merge(state.covariates,crosswalk[,1:2])
weights.for.post<-merge(weights.for.post,crosswalk[,1:2])


##============ready for MRP
library(arm)
##add state covariates to post-strat data
weights.for.post$race.gender<-factor(weights.for.post$race.gender,
                                     levels=c(1,2,3,4,5,6),
                                     labels=c("White.Male","White.Female",
                                              "Black.Male","Black.Female",
                                              "Other.Male","Other.Female"))
weights.for.post$age<-factor(weights.for.post$age,
                             levels=c(1,2,3,4,5),
                             labels=c("<18",
                                      "18 - 29",
                                      "30 - 44",
                                      "45 - 64",
                                      "65+"))
weights.for.post$educ<-factor(weights.for.post$educ,
                              levels=c(1,2,3,4),
                              labels=c("Less than HS",
                                       "HS",
                                       "Some college",
                                       "College grad"))


###percents total and percents for demographics
weights.for.post$percents<-weights.for.post$counts/weights.for.post$resp.numb




###start with death_penalty
death_penalty<-merge(death_penalty,state.covariates[,c(1,2,4,6,7)])

##create interactions--note that race.gender already is an interaction
death_penalty$race.gender.x.age<-interaction(death_penalty$race.gender,death_penalty$age)
death_penalty$race.gender.x.state<-interaction(death_penalty$race.gender,death_penalty$State)
death_penalty$race.gender.x.region<-interaction(death_penalty$race.gender,death_penalty$region)
death_penalty$age.x.state<-interaction(death_penalty$age,death_penalty$State)
death_penalty$age.x.region<-interaction(death_penalty$age,death_penalty$region)
death_penalty$state.x.year<-interaction(death_penalty$State,death_penalty$year)
death_penalty$race.gender.x.year<-interaction(death_penalty$race.gender,death_penalty$year)
death_penalty$age.x.year<-interaction(death_penalty$age,death_penalty$year)

#===run the MLM

library(lme4)
death_penalty$year.sq<-death_penalty$year*death_penalty$year
death_penalty$death<-as.numeric(death_penalty$death)

death_pen.mlm<-glmer(death~(1|race.gender)+
                       (1|race.gender.x.age)+
                       (1|race.gender.x.state)+
                       (1|race.gender.x.region)+
                       (1|race.gender.x.year)+
                       (1|age.x.year)+
                       (1|age.x.state)+
                       (1|age.x.region)+
                       (1|state.x.year)+
                       (1|educ)+
                       (1|age)+
                       black.percent+
                       pct.repub+
                       year+
                       year.sq,
                     data=death_penalty,
                     family=binomial(link="logit"))



##predict from weights.for.post.strat
##create interaciton ids in post.strat data
weights.for.post<-merge(weights.for.post,state.covariates[,c(1,2,4,6)],all.x=T)
weights.for.post$race.gender.x.age<-interaction(weights.for.post$race.gender,weights.for.post$age)
weights.for.post$race.gender.x.state<-interaction(weights.for.post$race.gender,weights.for.post$State)
weights.for.post$race.gender.x.region<-interaction(weights.for.post$race.gender,weights.for.post$region)
weights.for.post$race.gender.x.year<-interaction(weights.for.post$race.gender,weights.for.post$year)
weights.for.post$age.x.state<-interaction(weights.for.post$age,weights.for.post$State)
weights.for.post$age.x.region<-interaction(weights.for.post$age,weights.for.post$region)
weights.for.post$age.x.year<-interaction(weights.for.post$age,weights.for.post$year)
weights.for.post$state.x.year<-interaction(weights.for.post$State,weights.for.post$year)


##create wtihin group cell counts for disaggregation
weights.for.post$race<-"White"
weights.for.post$race[weights.for.post$race.gender=="Black.Male" | weights.for.post$race.gender=="Black.Female"]<-"Black"
weights.for.post$race[weights.for.post$race.gender=="Other.Male" | weights.for.post$race.gender=="Other.Female"]<-"Other"

#####------Start wtih raw totals


race.gender.eff<-ranef(death_pen.mlm)$race.gender[weights.for.post$race.gender,1]
race.gender.eff[which(is.na(race.gender.eff))]<-0 ##effects can be missing for states that were unsampled. Replace missings with 0 (e.g., w/ grand mean)

race.gender.age.eff<-ranef(death_pen.mlm)$race.gender.x.age[weights.for.post$race.gender.x.age,1]
race.gender.age.eff[which(is.na(race.gender.age.eff))]<-0

race.gender.year.eff<-ranef(death_pen.mlm)$race.gender.x.year[weights.for.post$race.gender.x.year,1]
race.gender.year.eff[which(is.na(race.gender.year.eff))]<-0


race.gender.state.eff<-ranef(death_pen.mlm)$race.gender.x.state[weights.for.post$race.gender.x.state,1]
race.gender.state.eff[which(is.na(race.gender.state.eff))]<-0

race.gender.region.eff<-ranef(death_pen.mlm)$race.gender.x.region[weights.for.post$race.gender.x.region,1]
race.gender.region.eff[which(is.na(race.gender.region.eff))]<-0

age.eff<-ranef(death_pen.mlm)$age[weights.for.post$age,1]
age.eff[which(is.na(age.eff))]<-0

age.state.eff<-ranef(death_pen.mlm)$age.x.state[weights.for.post$age.x.state,1]
age.state.eff[which(is.na(age.state.eff))]<-0


age.year.eff<-ranef(death_pen.mlm)$age.x.year[weights.for.post$age.x.year,1]
age.year.eff[which(is.na(age.year.eff))]<-0


age.region.eff<-ranef(death_pen.mlm)$age.x.region[weights.for.post$age.x.region,1]
age.region.eff[which(is.na(age.region.eff))]<-0

state.year.eff<-ranef(death_pen.mlm)$state.x.year[weights.for.post$state.x.year,1]
state.year.eff[which(is.na(state.year.eff))]<-0


educ.eff<-ranef(death_pen.mlm)$educ[weights.for.post$educ,1]
educ.eff[which(is.na(educ.eff))]<-0

bp.eff<-(fixef(death_pen.mlm)["black.percent"]*(weights.for.post$black.percent))
repub.eff<-(fixef(death_pen.mlm)["pct.repub"]*(weights.for.post$pct.repub))
year.eff<-(fixef(death_pen.mlm)["year"]*(weights.for.post$year))
yearsq.eff<-(fixef(death_pen.mlm)["year.sq"]*(weights.for.post$year.sq))


death_pen.pred<-arm::invlogit(fixef(death_pen.mlm)["(Intercept)"]+
                                race.gender.eff+
                                race.gender.age.eff+
                                race.gender.state.eff+
                                race.gender.region.eff+
                                race.gender.age.eff+
                                state.year.eff+
                                age.year.eff+
                                age.state.eff+
                                age.region.eff+
                                age.eff+
                                educ.eff+
                                bp.eff+
                                repub.eff+
                                year.eff+
                                yearsq.eff
)





###create weights 

###weight data
death_pen.weights<-death_pen.pred*weights.for.post$percents

#aggregate separately
weights.for.post$preds<-death_pen.weights
state.totals<-aggregate(preds~fipsst+year,weights,FUN=sum)



###Disaggregate by groups
#then post-stratify within groups

###create percents for race
black.counts.for.post<-weights.for.post[weights.for.post$race=="Black",]
black.counts.for.post<-aggregate(counts~fipsst+year+race,data=black.counts.for.post,FUN=sum)
colnames(black.counts.for.post)[4]<-"race.resp.numb"

white.counts.for.post<-weights.for.post[weights.for.post$race=="White",]
white.counts.for.post<-aggregate(counts~fipsst+year+race,data=white.counts.for.post,FUN=sum)
colnames(white.counts.for.post)[4]<-"race.resp.numb"

Other.counts.for.post<-weights.for.post[weights.for.post$race=="Other",]
Other.counts.for.post<-aggregate(counts~fipsst+year+race,data=Other.counts.for.post,FUN=sum)
colnames(Other.counts.for.post)[4]<-"race.resp.numb"


race.countts.for.post<-rbind(black.counts.for.post,white.counts.for.post,Other.counts.for.post)
weights.for.post<-merge(weights.for.post,race.countts.for.post)
weights.for.post$race.percents<-weights.for.post$counts/weights.for.post$race.resp.numb


##aggregate by race
weights.for.post$weights.race.weights<-weights.for.post$preds*weights.for.post$race.percents
race.estimates<-aggregate(weights.race.weights~year+fipsst+race,data=weights.for.post,FUN=sum)



##outcome disaggregated by race and overall
death_pen.race.totals<-race.estimates
death_pen.state.totals<-state.totals






###Dyad ratio algorithm implemented in standalone Wcalc5 software available from Stimson's website. 







###########################
# Primary analysis with ECM 
###########################




load("State control variables.RData")
load("Punitiveness by demographic 1972 - 2015.RData")
load("State crosswalk.RData")
load("Prison population 1970 - 2015.RData")

main.data<-state.data[,c(1:4,8:ncol(state.data))]
punitiveness.by.group$fipsst<-as.character(punitiveness.by.group$fipsst)
punitiveness.by.group[,3:ncol(punitiveness.by.group)]<-punitiveness.by.group[,3:ncol(punitiveness.by.group)]*100
main.data<-merge(main.data,punitiveness.by.group,all.x=T)
main.data<-main.data[main.data$year>=1969,]
prison.dat<-merge(main.data,prisoners[,-c(1)])


library(Amelia)##impute before merge---do this with 5 imputations when running final models
imputations<-amelia(prison.dat[,-c(2)],m=10,cs=c("State"),ts="year")


##function to convert imputation data to ECM data structure
ECM_process<-function(prison.dat){
  
  lag.dat<-prison.dat
  lag.dat$year<-lag.dat$year+1
  colnames(lag.dat)[3:ncol(lag.dat)]<-paste("lag.",colnames(prison.dat[3:ncol(lag.dat)]),sep="")
  
  
  diff.dat<-prison.dat
  colnames(diff.dat)[3:ncol(diff.dat)]<-paste("diff.",colnames(prison.dat)[3:ncol(diff.dat)],sep="")
  
  ecm.dat<-merge(diff.dat,lag.dat)
  ecm.dat$diff.pct.black<-ecm.dat$diff.pct.black-ecm.dat$lag.pct.black
  ecm.dat$diff.mood<-ecm.dat$diff.mood-ecm.dat$lag.mood
  ecm.dat$diff.Gini<-ecm.dat$diff.Gini-ecm.dat$lag.Gini
  ecm.dat$diff.prct.repub<-ecm.dat$diff.prct.repub-ecm.dat$lag.prct.repub
  ecm.dat$diff.unemp.rate<-ecm.dat$diff.unemp.rate-ecm.dat$lag.unemp.rate
  ecm.dat$diff.viol.crime.rate<-ecm.dat$diff.viol.crime.rate-ecm.dat$lag.viol.crime.rate
  ecm.dat$diff.punitiveness.blacks<-ecm.dat$diff.punitiveness.blacks-ecm.dat$lag.punitiveness.blacks
  ecm.dat$diff.punitiveness.whites<-ecm.dat$diff.punitiveness.whites-ecm.dat$lag.punitiveness.whites
  ecm.dat$diff.punitiveness.Males<-ecm.dat$diff.punitiveness.Males-ecm.dat$lag.punitiveness.Males
  ecm.dat$diff.punitiveness.Females<-ecm.dat$diff.punitiveness.Females-ecm.dat$lag.punitiveness.Females
  ecm.dat$`diff.punitiveness.65+`<-ecm.dat$`diff.punitiveness.65+`-ecm.dat$`lag.punitiveness.65+`
  ecm.dat$diff.punitiveness.45<-ecm.dat$diff.punitiveness.45-ecm.dat$lag.punitiveness.45
  ecm.dat$diff.punitiveness.18<-ecm.dat$diff.punitiveness.18-ecm.dat$lag.punitiveness.18
  ecm.dat$diff.punitiveness.30<-ecm.dat$diff.punitiveness.30-ecm.dat$lag.punitiveness.30
  ecm.dat$diff.incar.rate<-ecm.dat$diff.incar.rate-ecm.dat$lag.incar.rate
  
  return(ecm.dat)
  
}


##loop over imputations
for(i in 1:10){
imp.dat<-imputations$imputations[[i]]


##rduce to only variables I wil luse
prison.dat<-imp.dat[,-c(3,4,6,7,21,22)]


ecm.dat<-ECM_process(prison.dat)


library(lmtest)
library(sandwich)
library(plm)
tdat<-pdata.frame(ecm.dat,index=c("State","year"))
tdat<-tdat[,-c(1:2)]

for(t in 1:ncol(tdat)){
  print(colnames(tdat)[t])
 print(purtest(tdat[,t]))
}

#stationry variables:
#lag.mood
#lag.unemp.rate
#lag.pct.black
#all punitiveness variables

out<-purtest(tdat[,-c(1:2)],test="levinlin") ##yes, stationary model
                                           
                                                      



###run model
model1<-lm(diff.incar.rate~
             #error correction rate
             lag.incar.rate+
             #transient terms
             #diff.pct.black+
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.blacks+
             lag.punitiveness.whites,
           data=ecm.dat)


coef_output<-coeftest(model1, vcov = vcovHC(model1, type="HC1"))

library(car)
car::vif(model1) #check collinearity


###store coefficients
if(i==1){
  coef.mat<-coef_output[,1]
  se.mat<-coef_output[,2]
}else{
  coef.mat2<-coef_output[,1]
  se.mat2<-coef_output[,2]
  
  coef.mat<-rbind(coef.mat,coef.mat2)
  se.mat<-rbind(se.mat,se.mat2)
}

    if(i==1){
          LRM_coef<-cbind(deltaMethod(model1,"lag.pct.black/1-lag.incar.rate")[,1],
        deltaMethod(model1,"lag.mood/1-lag.incar.rate")[,1],
        deltaMethod(model1,"lag.Gini/1-lag.incar.rate")[,1],
        deltaMethod(model1,"lag.prct.repub/1-lag.incar.rate")[,1],
        deltaMethod(model1,"lag.unemp.rate/1-lag.incar.rate")[,1],
        deltaMethod(model1,"lag.viol.crime.rate/1-lag.incar.rate")[,1],
        deltaMethod(model1,"lag.punitiveness.blacks/1-lag.incar.rate")[,1],
        deltaMethod(model1,"lag.punitiveness.whites/1-lag.incar.rate")[,1])
          
        LRM_SE<-cbind(deltaMethod(model1,"lag.pct.black/1-lag.incar.rate")[,2],
                          deltaMethod(model1,"lag.mood/1-lag.incar.rate")[,2],
                          deltaMethod(model1,"lag.Gini/1-lag.incar.rate")[,2],
                          deltaMethod(model1,"lag.prct.repub/1-lag.incar.rate")[,2],
                          deltaMethod(model1,"lag.unemp.rate/1-lag.incar.rate")[,2],
                          deltaMethod(model1,"lag.viol.crime.rate/1-lag.incar.rate")[,2],
                          deltaMethod(model1,"lag.punitiveness.blacks/1-lag.incar.rate")[,2],
                          deltaMethod(model1,"lag.punitiveness.whites/1-lag.incar.rate")[,2])
    }else{
      
      LRM_coef2<-cbind(deltaMethod(model1,"lag.pct.black/1-lag.incar.rate")[,1],
                      deltaMethod(model1,"lag.mood/1-lag.incar.rate")[,1],
                      deltaMethod(model1,"lag.Gini/1-lag.incar.rate")[,1],
                      deltaMethod(model1,"lag.prct.repub/1-lag.incar.rate")[,1],
                      deltaMethod(model1,"lag.unemp.rate/1-lag.incar.rate")[,1],
                      deltaMethod(model1,"lag.viol.crime.rate/1-lag.incar.rate")[,1],
                      deltaMethod(model1,"lag.punitiveness.blacks/1-lag.incar.rate")[,1],
                      deltaMethod(model1,"lag.punitiveness.whites/1-lag.incar.rate")[,1])
      
      LRM_SE2<-cbind(deltaMethod(model1,"lag.pct.black/1-lag.incar.rate")[,2],
                    deltaMethod(model1,"lag.mood/1-lag.incar.rate")[,2],
                    deltaMethod(model1,"lag.Gini/1-lag.incar.rate")[,2],
                    deltaMethod(model1,"lag.prct.repub/1-lag.incar.rate")[,2],
                    deltaMethod(model1,"lag.unemp.rate/1-lag.incar.rate")[,2],
                    deltaMethod(model1,"lag.viol.crime.rate/1-lag.incar.rate")[,2],
                    deltaMethod(model1,"lag.punitiveness.blacks/1-lag.incar.rate")[,2],
                    deltaMethod(model1,"lag.punitiveness.whites/1-lag.incar.rate")[,2])
      
      
      LRM_coef<-rbind(LRM_coef,LRM_coef2)
      LRM_SE<-rbind(LRM_SE,LRM_SE2)
      
      
    }


}
pooled_coefs1<-mi.meld(coef.mat,se.mat)
pooled_LRM1<-mi.meld(LRM_coef,LRM_SE)



###model 2


for(i in 1:10){
  imp.dat<-imputations$imputations[[i]]
  
  
  ##rduce to only variables I wil luse
  prison.dat<-imp.dat[,-c(3,4,6,7,21,22)]
  
  
  ecm.dat<-ECM_process(prison.dat)

###run model
model2<-lm(diff.incar.rate~
             #error correction rate
             lag.incar.rate+
             #transient terms
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.Males+
             lag.punitiveness.Females,
           data=ecm.dat)



coef_output<-coeftest(model2, vcov = vcovHC(model2, type="HC1"))

library(car)
car::vif(model2) #check collin


###store coefficients
if(i==1){
  coef.mat<-coef_output[,1]
  se.mat<-coef_output[,2]
}else{
  coef.mat2<-coef_output[,1]
  se.mat2<-coef_output[,2]
  
  coef.mat<-rbind(coef.mat,coef.mat2)
  se.mat<-rbind(se.mat,se.mat2)
}

if(i==1){
  LRM_coef<-cbind(deltaMethod(model2,"lag.pct.black/1-lag.incar.rate")[,1],
                  deltaMethod(model2,"lag.mood/1-lag.incar.rate")[,1],
                  deltaMethod(model2,"lag.Gini/1-lag.incar.rate")[,1],
                  deltaMethod(model2,"lag.prct.repub/1-lag.incar.rate")[,1],
                  deltaMethod(model2,"lag.unemp.rate/1-lag.incar.rate")[,1],
                  deltaMethod(model2,"lag.viol.crime.rate/1-lag.incar.rate")[,1],
                  deltaMethod(model2,"lag.punitiveness.Males/1-lag.incar.rate")[,1],
                  deltaMethod(model2,"lag.punitiveness.Females/1-lag.incar.rate")[,1])
  
  LRM_SE<-cbind(deltaMethod(model2,"lag.pct.black/1-lag.incar.rate")[,2],
                deltaMethod(model2,"lag.mood/1-lag.incar.rate")[,2],
                deltaMethod(model2,"lag.Gini/1-lag.incar.rate")[,2],
                deltaMethod(model2,"lag.prct.repub/1-lag.incar.rate")[,2],
                deltaMethod(model2,"lag.unemp.rate/1-lag.incar.rate")[,2],
                deltaMethod(model2,"lag.viol.crime.rate/1-lag.incar.rate")[,2],
                deltaMethod(model2,"lag.punitiveness.Males/1-lag.incar.rate")[,2],
                deltaMethod(model2,"lag.punitiveness.Females/1-lag.incar.rate")[,2])
}else{
  
  LRM_coef2<-cbind(deltaMethod(model2,"lag.pct.black/1-lag.incar.rate")[,1],
                   deltaMethod(model2,"lag.mood/1-lag.incar.rate")[,1],
                   deltaMethod(model2,"lag.Gini/1-lag.incar.rate")[,1],
                   deltaMethod(model2,"lag.prct.repub/1-lag.incar.rate")[,1],
                   deltaMethod(model2,"lag.unemp.rate/1-lag.incar.rate")[,1],
                   deltaMethod(model2,"lag.viol.crime.rate/1-lag.incar.rate")[,1],
                   deltaMethod(model2,"lag.punitiveness.Males/1-lag.incar.rate")[,1],
                   deltaMethod(model2,"lag.punitiveness.Females/1-lag.incar.rate")[,1])
  
  LRM_SE2<-cbind(deltaMethod(model2,"lag.pct.black/1-lag.incar.rate")[,2],
                 deltaMethod(model2,"lag.mood/1-lag.incar.rate")[,2],
                 deltaMethod(model2,"lag.Gini/1-lag.incar.rate")[,2],
                 deltaMethod(model2,"lag.prct.repub/1-lag.incar.rate")[,2],
                 deltaMethod(model2,"lag.unemp.rate/1-lag.incar.rate")[,2],
                 deltaMethod(model2,"lag.viol.crime.rate/1-lag.incar.rate")[,2],
                 deltaMethod(model2,"lag.punitiveness.Males/1-lag.incar.rate")[,2],
                 deltaMethod(model2,"lag.punitiveness.Females/1-lag.incar.rate")[,2])
  
  
  LRM_coef<-rbind(LRM_coef,LRM_coef2)
  LRM_SE<-rbind(LRM_SE,LRM_SE2)
  
  
}




}

pooled_coefs2<-mi.meld(coef.mat,se.mat)
pooled_LRM2<-mi.meld(LRM_coef,LRM_SE)




ecm.dat$lag.punitiveness.65<-ecm.dat$`lag.punitiveness.65+`


#model 3


for(i in 1:10){
  imp.dat<-imputations$imputations[[i]]
  
  ##rduce to only variables I wil luse
  prison.dat<-imp.dat[,-c(3,4,6,7,21,22)]
  
  
  ecm.dat<-ECM_process(prison.dat)
  
###run model
model3<-lm(diff.incar.rate~
             #error correction rate
             lag.incar.rate+
             #transient terms
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.18+
             lag.punitiveness.30+
             lag.punitiveness.45+
             `lag.punitiveness.65+`,
           data=ecm.dat)




coef_output<-coeftest(model3, vcov = vcovHC(model3, type="HC1"))

library(car)
car::vif(model3) #check collin


###store coefficients
if(i==1){
  coef.mat<-coef_output[,1]
  se.mat<-coef_output[,2]
}else{
  coef.mat2<-coef_output[,1]
  se.mat2<-coef_output[,2]
  
  coef.mat<-rbind(coef.mat,coef.mat2)
  se.mat<-rbind(se.mat,se.mat2)
}

if(i==1){
  LRM_coef<-cbind(deltaMethod(model3,"lag.pct.black/1-lag.incar.rate")[,1],
                  deltaMethod(model3,"lag.mood/1-lag.incar.rate")[,1],
                  deltaMethod(model3,"lag.Gini/1-lag.incar.rate")[,1],
                  deltaMethod(model3,"lag.prct.repub/1-lag.incar.rate")[,1],
                  deltaMethod(model3,"lag.unemp.rate/1-lag.incar.rate")[,1],
                  deltaMethod(model3,"lag.viol.crime.rate/1-lag.incar.rate")[,1],
                  deltaMethod(model3,"lag.punitiveness.18/1-lag.incar.rate")[,1],
                  deltaMethod(model3,"lag.punitiveness.30/1-lag.incar.rate")[,1],
                  deltaMethod(model3,"lag.punitiveness.45/1-lag.incar.rate")[,1],
                  deltaMethod(model3,"lag.punitiveness.65/1-lag.incar.rate")[,1])
                  
  LRM_SE<-cbind(deltaMethod(model3,"lag.pct.black/1-lag.incar.rate")[,2],
                deltaMethod(model3,"lag.mood/1-lag.incar.rate")[,2],
                deltaMethod(model3,"lag.Gini/1-lag.incar.rate")[,2],
                deltaMethod(model3,"lag.prct.repub/1-lag.incar.rate")[,2],
                deltaMethod(model3,"lag.unemp.rate/1-lag.incar.rate")[,2],
                deltaMethod(model3,"lag.viol.crime.rate/1-lag.incar.rate")[,2],
                deltaMethod(model3,"lag.punitiveness.18/1-lag.incar.rate")[,2],
                deltaMethod(model3,"lag.punitiveness.30/1-lag.incar.rate")[,2],
                deltaMethod(model3,"lag.punitiveness.45/1-lag.incar.rate")[,2],
                deltaMethod(model3,"lag.punitiveness.65/1-lag.incar.rate")[,2])
}else{
  
  LRM_coef2<-cbind(deltaMethod(model3,"lag.pct.black/1-lag.incar.rate")[,1],
                   deltaMethod(model3,"lag.mood/1-lag.incar.rate")[,1],
                   deltaMethod(model3,"lag.Gini/1-lag.incar.rate")[,1],
                   deltaMethod(model3,"lag.prct.repub/1-lag.incar.rate")[,1],
                   deltaMethod(model3,"lag.unemp.rate/1-lag.incar.rate")[,1],
                   deltaMethod(model3,"lag.viol.crime.rate/1-lag.incar.rate")[,1],
                   deltaMethod(model3,"lag.punitiveness.18/1-lag.incar.rate")[,1],
                   deltaMethod(model3,"lag.punitiveness.30/1-lag.incar.rate")[,1],
                   deltaMethod(model3,"lag.punitiveness.45/1-lag.incar.rate")[,1],
                   deltaMethod(model3,"lag.punitiveness.65/1-lag.incar.rate")[,1])
  
  LRM_SE2<-cbind(deltaMethod(model3,"lag.pct.black/1-lag.incar.rate")[,2],
                 deltaMethod(model3,"lag.mood/1-lag.incar.rate")[,2],
                 deltaMethod(model3,"lag.Gini/1-lag.incar.rate")[,2],
                 deltaMethod(model3,"lag.prct.repub/1-lag.incar.rate")[,2],
                 deltaMethod(model3,"lag.unemp.rate/1-lag.incar.rate")[,2],
                 deltaMethod(model3,"lag.viol.crime.rate/1-lag.incar.rate")[,2],
                 deltaMethod(model3,"lag.punitiveness.18/1-lag.incar.rate")[,2],
                 deltaMethod(model3,"lag.punitiveness.30/1-lag.incar.rate")[,2],
                 deltaMethod(model3,"lag.punitiveness.45/1-lag.incar.rate")[,2],
                 deltaMethod(model3,"lag.punitiveness.65/1-lag.incar.rate")[,2])
  
  
  LRM_coef<-rbind(LRM_coef,LRM_coef2)
  LRM_SE<-rbind(LRM_SE,LRM_SE2)
  
  
}




}

pooled_coefs3<-mi.meld(coef.mat,se.mat)
pooled_LRM3<-mi.meld(LRM_coef,LRM_SE)






##now with ridge regression
library(ridge)

library(lmridge)

##model 4


for(i in 1:10){
  imp.dat<-imputations$imputations[[i]]
  
  ##rduce to only variables I wil luse
  prison.dat<-imp.dat[,-c(3,4,6,7,21,22)]
  
  
  ecm.dat<-ECM_process(prison.dat)
  
###run model
model4<-linearRidge(diff.incar.rate~
                  #error correction rate
                  lag.incar.rate+
                  #transient terms
                  #diff.pct.black+
                  diff.Gini+
                  diff.prct.repub+
                  diff.viol.crime.rate+
                  #equlibrium terms
                  lag.pct.black+
                  lag.mood+
                  lag.Gini+
                  lag.prct.repub+
                  lag.unemp.rate+
                  lag.viol.crime.rate+
                  lag.punitiveness.blacks+
                  lag.punitiveness.whites,
                data=ecm.dat,scaling="none",
                lambda=.344) ##lambda pre-estimated using the "lmridge" function and fixed here to decrease estimation times


###store coefficients
if(i==1){
  coef.mat<-t(as.matrix(coef(model4)))
  se.mat<-t(as.matrix(sqrt(diag(vcov(model4)))))
}else{
  coef.mat2<-t(as.matrix(coef(model4)))
  se.mat2<-t(as.matrix(sqrt(diag(vcov(model4)))))
  
  coef.mat<-rbind(coef.mat,coef.mat2)
  se.mat<-rbind(se.mat,se.mat2)
}

##again for LRM

model4lrm<-lmridge(diff.incar.rate~
             #error correction rate
             lag.incar.rate+
             #transient terms
             #diff.pct.black+
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.blacks+
             lag.punitiveness.whites,
           data=ecm.dat,
           model4$lambda[1])

summ<-summary(model4lrm)
varcov1<-vcov(model4lrm)[[1]]

if(i==1){
  LRM_coef<-cbind(deltaMethod(model4lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov1)[,1],
                  deltaMethod(model4lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov1)[,1],
                  deltaMethod(model4lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov1)[,1],
                  deltaMethod(model4lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov1)[,1],
                  deltaMethod(model4lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov1)[,1],
                  deltaMethod(model4lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov1)[,1],
                  deltaMethod(model4lrm$coef[,1],"lag.punitiveness.blacks/1-lag.incar.rate",vcov=varcov1)[,1],
                  deltaMethod(model4lrm$coef[,1],"lag.punitiveness.whites/1-lag.incar.rate",vcov=varcov1)[,1])
  
  LRM_SE<-cbind(deltaMethod(model4lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov1)[,2],
                deltaMethod(model4lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov1)[,2],
                deltaMethod(model4lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov1)[,2],
                deltaMethod(model4lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov1)[,2],
                deltaMethod(model4lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov1)[,2],
                deltaMethod(model4lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov1)[,2],
                deltaMethod(model4lrm$coef[,1],"lag.punitiveness.blacks/1-lag.incar.rate",vcov=varcov1)[,2],
                deltaMethod(model4lrm$coef[,1],"lag.punitiveness.whites/1-lag.incar.rate",vcov=varcov1)[,2])
}else{
  
  LRM_coef2<-cbind(deltaMethod(model4lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov1)[,1],
                   deltaMethod(model4lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov1)[,1],
                   deltaMethod(model4lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov1)[,1],
                   deltaMethod(model4lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov1)[,1],
                   deltaMethod(model4lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov1)[,1],
                   deltaMethod(model4lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov1)[,1],
                   deltaMethod(model4lrm$coef[,1],"lag.punitiveness.blacks/1-lag.incar.rate",vcov=varcov1)[,1],
                   deltaMethod(model4lrm$coef[,1],"lag.punitiveness.whites/1-lag.incar.rate",vcov=varcov1)[,1])
  
  LRM_SE2<-cbind(deltaMethod(model4lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov1)[,2],
                 deltaMethod(model4lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov1)[,2],
                 deltaMethod(model4lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov1)[,2],
                 deltaMethod(model4lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov1)[,2],
                 deltaMethod(model4lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov1)[,2],
                 deltaMethod(model4lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov1)[,2],
                 deltaMethod(model4lrm$coef[,1],"lag.punitiveness.blacks/1-lag.incar.rate",vcov=varcov1)[,2],
                 deltaMethod(model4lrm$coef[,1],"lag.punitiveness.whites/1-lag.incar.rate",vcov=varcov1)[,2])
  
  
  LRM_coef<-rbind(LRM_coef,LRM_coef2)
  LRM_SE<-rbind(LRM_SE,LRM_SE2)
  
  
}




}



pooled_coefs4<-mi.meld(coef.mat,se.mat)
pooled_LRM4<-mi.meld(LRM_coef,LRM_SE)


###model 5


for(i in 1:10){
  imp.dat<-imputations$imputations[[i]]
  
  ##rduce to only variables I wil luse
  prison.dat<-imp.dat[,-c(3,4,6,7,21,22)]
  
  
  ecm.dat<-ECM_process(prison.dat)
###run model
model5<-linearRidge(diff.incar.rate~
                      #error correction rate
                      lag.incar.rate+
                      #transient terms
                      #diff.pct.black+
                      diff.Gini+
                      diff.prct.repub+
                      diff.viol.crime.rate+
                      #equlibrium terms
                      lag.pct.black+
                      lag.mood+
                      lag.Gini+
                      lag.prct.repub+
                      lag.unemp.rate+
                      lag.viol.crime.rate+
                      lag.punitiveness.Males+
                      lag.punitiveness.Females,
                    data=ecm.dat,scaling="none",lambda=.244)


###store coefficients
if(i==1){
  coef.mat<-t(as.matrix(coef(model5)))
  se.mat<-t(as.matrix(sqrt(diag(vcov(model5)))))
}else{
  coef.mat2<-t(as.matrix(coef(model5)))
  se.mat2<-t(as.matrix(sqrt(diag(vcov(model5)))))
  
  coef.mat<-rbind(coef.mat,coef.mat2)
  se.mat<-rbind(se.mat,se.mat2)
}


##again for LRM

model5lrm<-lmridge(diff.incar.rate~
                     #error correction rate
                     lag.incar.rate+
                     #transient terms
                     #diff.pct.black+
                     diff.Gini+
                     diff.prct.repub+
                     diff.viol.crime.rate+
                     #equlibrium terms
                     lag.pct.black+
                     lag.mood+
                     lag.Gini+
                     lag.prct.repub+
                     lag.unemp.rate+
                     lag.viol.crime.rate+
                     lag.punitiveness.Males+
                     lag.punitiveness.Females,
                   data=ecm.dat,
                   K=model5$lambda[1])

summ2<-summary(model5lrm)

varcov2<-vcov(model5lrm)[[1]]
if(i==1){
  LRM_coef<-cbind(deltaMethod(model5lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model5lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model5lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model5lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model5lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model5lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model5lrm$coef[,1],"lag.punitiveness.Males/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model5lrm$coef[,1],"lag.punitiveness.Females/1-lag.incar.rate",vcov=varcov2)[,1])
  
  LRM_SE<-cbind(deltaMethod(model5lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model5lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model5lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model5lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model5lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model5lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model5lrm$coef[,1],"lag.punitiveness.Males/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model5lrm$coef[,1],"lag.punitiveness.Females/1-lag.incar.rate",vcov=varcov2)[,2])
}else{
  
  LRM_coef2<-cbind(deltaMethod(model5lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model5lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model5lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model5lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model5lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model5lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model5lrm$coef[,1],"lag.punitiveness.Males/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model5lrm$coef[,1],"lag.punitiveness.Females/1-lag.incar.rate",vcov=varcov2)[,1])
  
  LRM_SE2<-cbind(deltaMethod(model5lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model5lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model5lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model5lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model5lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model5lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model5lrm$coef[,1],"lag.punitiveness.Males/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model5lrm$coef[,1],"lag.punitiveness.Females/1-lag.incar.rate",vcov=varcov2)[,2])
  
  
  LRM_coef<-rbind(LRM_coef,LRM_coef2)
  LRM_SE<-rbind(LRM_SE,LRM_SE2)
  
  
}




}


pooled_coefs5<-mi.meld(coef.mat,se.mat)
pooled_LRM5<-mi.meld(LRM_coef,LRM_SE)



###model 6


for(i in 1:10){
  imp.dat<-imputations$imputations[[i]]
  
  ##rduce to only variables I wil luse
  prison.dat<-imp.dat[,-c(3,4,6,7,21,22)]
  
  
  ecm.dat<-ECM_process(prison.dat)

###run model
model6<-linearRidge(diff.incar.rate~
                      #error correction rate
                      lag.incar.rate+
                      #transient terms
                      #diff.pct.black+
                      diff.Gini+
                      diff.prct.repub+
                      diff.viol.crime.rate+
                      #equlibrium terms
                      lag.pct.black+
                      lag.mood+
                      lag.Gini+
                      lag.prct.repub+
                      lag.unemp.rate+
                      lag.viol.crime.rate+
                      lag.punitiveness.18+
                      lag.punitiveness.30+
                      lag.punitiveness.45+
                      lag.punitiveness.65,
                    data=ecm.dat,scaling="none",lambda=.343)



##again for LRM

model6lrm<-lmridge(diff.incar.rate~
                     #error correction rate
                     lag.incar.rate+
                     #transient terms
                     #diff.pct.black+
                     diff.Gini+
                     diff.prct.repub+
                     diff.viol.crime.rate+
                     #equlibrium terms
                     lag.pct.black+
                     lag.mood+
                     lag.Gini+
                     lag.prct.repub+
                     lag.unemp.rate+
                     lag.viol.crime.rate+
                     lag.punitiveness.18+
                     lag.punitiveness.30+
                     lag.punitiveness.45+
                     lag.punitiveness.65,
                   data=ecm.dat,
                   K=model6$lambda[1])

summ3<-summary(model6lrm)

varcov2<-vcov(model6lrm)[[1]]
if(i==1){
  LRM_coef<-cbind(deltaMethod(model6lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model6lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model6lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model6lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model6lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model6lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model6lrm$coef[,1],"lag.punitiveness.18/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model6lrm$coef[,1],"lag.punitiveness.30/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model6lrm$coef[,1],"lag.punitiveness.45/1-lag.incar.rate",vcov=varcov2)[,1],
                  deltaMethod(model6lrm$coef[,1],"lag.punitiveness.65/1-lag.incar.rate",vcov=varcov2)[,1])
                  
  LRM_SE<-cbind(deltaMethod(model6lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model6lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model6lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model6lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model6lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model6lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model6lrm$coef[,1],"lag.punitiveness.18/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model6lrm$coef[,1],"lag.punitiveness.30/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model6lrm$coef[,1],"lag.punitiveness.45/1-lag.incar.rate",vcov=varcov2)[,2],
                deltaMethod(model6lrm$coef[,1],"lag.punitiveness.65/1-lag.incar.rate",vcov=varcov2)[,2])
  
}else{
  
  LRM_coef2<-cbind(deltaMethod(model6lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model6lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model6lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model6lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model6lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model6lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model6lrm$coef[,1],"lag.punitiveness.18/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model6lrm$coef[,1],"lag.punitiveness.30/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model6lrm$coef[,1],"lag.punitiveness.45/1-lag.incar.rate",vcov=varcov2)[,1],
                   deltaMethod(model6lrm$coef[,1],"lag.punitiveness.65/1-lag.incar.rate",vcov=varcov2)[,1])
  
  
  LRM_SE2<-cbind(deltaMethod(model6lrm$coef[,1],"lag.pct.black/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model6lrm$coef[,1],"lag.mood/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model6lrm$coef[,1],"lag.Gini/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model6lrm$coef[,1],"lag.prct.repub/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model6lrm$coef[,1],"lag.unemp.rate/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model6lrm$coef[,1],"lag.viol.crime.rate/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model6lrm$coef[,1],"lag.punitiveness.18/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model6lrm$coef[,1],"lag.punitiveness.30/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model6lrm$coef[,1],"lag.punitiveness.45/1-lag.incar.rate",vcov=varcov2)[,2],
                 deltaMethod(model6lrm$coef[,1],"lag.punitiveness.65/1-lag.incar.rate",vcov=varcov2)[,2])
  
  
  
  LRM_coef<-rbind(LRM_coef,LRM_coef2)
  LRM_SE<-rbind(LRM_SE,LRM_SE2)
  
  
}




}

pooled_coefs6<-mi.meld(coef.mat,se.mat)
pooled_LRM6<-mi.meld(LRM_coef,LRM_SE)















#########################
#  decomposition using ECM predictions
###############


###aggregate to time period
decomp.dat1<-model1$model
decomp.dat1$year<-ecm.dat$year
decomp.dat1$State<-ecm.dat$State

model1_pred<-decomp.dat1
model1_pred$model1pred<-predict(model1)
model1_pred<-aggregate(model1pred~year,FUN=mean,data=model1_pred)

decomp.dat2<-model2$model
decomp.dat2$year<-ecm.dat$year
decomp.dat2$State<-ecm.dat$State

model2_pred<-decomp.dat2
model2_pred$model2pred<-predict(model2)
model2_pred<-aggregate(model2pred~year,FUN=mean,data=model2_pred)

decomp.dat3<-model3$model
decomp.dat3$year<-ecm.dat$year
decomp.dat3$State<-ecm.dat$State

model3_pred<-decomp.dat3
model3_pred$model3pred<-predict(model3)
model3_pred<-aggregate(model3pred~year,FUN=mean,data=model3_pred)



###fix values at 1970 values

pw.1970<-imp.dat[imp.dat$year==1970,c(2,14)]
colnames(pw.1970)[2]<-"white.pun.1970"
pw.1970[,2]<-pw.1970[,2]-5
wh.pun<-merge(decomp.dat1,pw.1970)
wh.pun<-wh.pun[,-c(1:2,14:15)]

wh.pun$pred <- apply(wh.pun, 1, function(x) t(x) %*% model1$coefficients[-c(1)])
wh.pun$pred<-wh.pun$pred+model1$coefficients[1]
wh.pun$year<-decomp.dat$year
wh.pun<-aggregate(pred~year,FUN=mean,data=wh.pun)


pw.1970<-imp.dat[imp.dat$year==1970,c(2,15)]
colnames(pw.1970)[2]<-"male.pun.1970"
pw.1970[,2]<-pw.1970[,2]-5

men.pun<-merge(decomp.dat2,pw.1970)
men.pun$lag.punitiveness.Males<-men.pun$male.pun.1970
men.pun<-men.pun[,-c(1:2,15,16)]

men.pun$pred <- apply(men.pun, 1, function(x) t(x) %*% model2$coefficients[-c(1)])
men.pun$pred<-men.pun$pred+model2$coefficients[1]
men.pun$year<-decomp.dat2$year
men.pun<-aggregate(pred~year,FUN=mean,data=men.pun)




pw.1970<-imp.dat[imp.dat$year==1970,c(2,20)]
colnames(pw.1970)[2]<-"middle.age.pun.1970"
pw.1970[,2]<-pw.1970[,2]-5

middle.age.pun<-merge(decomp.dat3,pw.1970)
middle.age.pun$lag.punitiveness.30<-middle.age.pun$middle.age.pun.1970
middle.age.pun<-middle.age.pun[,-c(1:2,17,18)]

middle.age.pun$pred <- apply(middle.age.pun, 1, function(x) t(x) %*% model3$coefficients[-c(1)])
middle.age.pun$pred<-middle.age.pun$pred+model3$coefficients[1]
middle.age.pun$year<-decomp.dat2$year
middle.age.pun<-aggregate(pred~year,FUN=mean,data=middle.age.pun)


##sum up change scores to provide predicted incarceration rates
library(dplyr)
wh.pun$pred<-cumsum(wh.pun$pred)
wh.pun$Model<-"Whites' punitiveness"
wh.pun$pred.type<-"Held at 1970"

gg.dat<-model1_pred
colnames(gg.dat)[2]<-"pred"
gg.dat$pred<-cumsum(gg.dat$pred)
gg.dat$Model<-"Whites' punitiveness"
gg.dat$pred.type<-"Predicted"

gg.dat<-rbind(gg.dat,wh.pun)


men.pun$pred<-cumsum(men.pun$pred)
men.pun$Model<-"Mens' punitiveness"
men.pun$pred.type<-"Held at 1970"

gg.dat2<-model2_pred
colnames(gg.dat2)[2]<-"pred"
gg.dat2$pred<-cumsum(gg.dat2$pred)
gg.dat2$Model<-"Mens' punitiveness"
gg.dat2$pred.type<-"Predicted"

gg.dat<-rbind(gg.dat,men.pun,gg.dat2)



middle.age.pun$pred<-cumsum(middle.age.pun$pred)
middle.age.pun$Model<-"30 - 44 year-\nold punitiveness"
middle.age.pun$pred.type<-"Held at 1970"

gg.dat2<-model3_pred
colnames(gg.dat2)[2]<-"pred"
gg.dat2$pred<-cumsum(gg.dat2$pred)
gg.dat2$Model<-"30 - 44 year-\nold punitiveness"
gg.dat2$pred.type<-"Predicted"

gg.dat<-rbind(gg.dat,middle.age.pun,gg.dat2)
gg.dat$pred<-gg.dat$pred+mean(imp.dat$incar.rate[imp.dat$year==1970])




gg.dat2<-aggregate(incar.rate~year,FUN=mean,data=prison.dat)
gg.dat2$Model<-"30 - 44 year-\nold punitiveness"
gg.dat2$pred.type<-"Observed"

gg.dat3<-gg.dat2
gg.dat2<-aggregate(incar.rate~year,FUN=mean,data=prison.dat)
gg.dat2$Model<-"Mens' punitiveness"
gg.dat2$pred.type<-"Observed"
gg.dat3<-rbind(gg.dat3,gg.dat2)

gg.dat2<-aggregate(incar.rate~year,FUN=mean,data=prison.dat)
gg.dat2$Model<-"Whites' punitiveness"
gg.dat2$pred.type<-"Observed"
gg.dat3<-rbind(gg.dat3,gg.dat2)

colnames(gg.dat3)[2]<-"pred"
gg.dat<-rbind(gg.dat3,gg.dat)




##ggplot
library(ggplot2)
jpeg("Decomposition.jpeg", width=10, height =6, units="in", res=1200)

ggplot(gg.dat,aes(x=year,y=pred,linetype=pred.type),size=1.2)+geom_line()+
  facet_wrap(.~Model)+ylab("Incarceration rate\nper 100,000")+
  theme_bw()+
  theme(legend.title=element_blank())+
  xlab("")

dev.off()



###do violent crime

pw.1970<-imp.dat[imp.dat$year==1970,c(2,12)]
colnames(pw.1970)[2]<-"vcrime"
wh.pun<-merge(decomp.dat1,pw.1970)
wh.pun$lag.viol.crime.rate<-wh.pun$vcrime
wh.pun<-wh.pun[,-c(1:2,15:16)]

wh.pun$pred <- apply(wh.pun, 1, function(x) t(x) %*% model1$coefficients[-c(1)])
wh.pun$pred<-wh.pun$pred+model1$coefficients[1]
wh.pun$year<-decomp.dat$year
wh.pun<-aggregate(pred~year,FUN=mean,data=wh.pun)
cumsum(wh.pun$pred)+mean(imp.dat$incar.rate[imp.dat$year==1970])

















####mediation analysis
library(mediation)
###run model



model7<-lm(diff.incar.rate~
                 year+
                 #error correction rate
                 lag.incar.rate+
                 #transient terms
                 #diff.pct.black+
                 diff.Gini+
                 diff.prct.repub+
                 diff.viol.crime.rate+
                 #equlibrium terms
                 lag.pct.black+
                 lag.mood+
                 lag.Gini+
                 lag.prct.repub+
                 lag.unemp.rate+
                 lag.viol.crime.rate+
                 lag.punitiveness.blacks,
               data=ecm.dat)

coeftest(model7, vcov = vcovHC(model7, type="HC1"))


model8<-lm(diff.incar.rate~
             year+
             #error correction rate
             lag.incar.rate+
             #transient terms
             #diff.pct.black+
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.whites,
           data=ecm.dat)

coeftest(model8, vcov = vcovHC(model8, type="HC1"))



model1.med<-lm(lag.punitiveness.whites~
             #error correction rate
             lag.incar.rate+
             #transient terms
             #diff.pct.black+
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.blacks,
           data=ecm.dat)

coeftest(model1.med, vcov = vcovHC(model1.med, type="HC"))

set.seed(21093)
whites.med<-mediate(model.m=model1.med,
                    model.y=model1,
                    treat="lag.punitiveness.blacks",
                    mediator="lag.punitiveness.whites",
                    long=FALSE,
                    robustSE=T,sims=10000)







model9<-lm(diff.incar.rate~
             year+
             #error correction rate
             lag.incar.rate+
             #transient terms
             #diff.pct.black+
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.Females,
           data=ecm.dat)

coeftest(model9, vcov = vcovHC(model9, type="HC1"))


model10<-lm(diff.incar.rate~
             year+
             #error correction rate
             lag.incar.rate+
             #transient terms
             #diff.pct.black+
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.Males,
           data=ecm.dat)

coeftest(model10, vcov = vcovHC(model10, type="HC1"))


###run model
model2.med<-lm(lag.punitiveness.Males~
             #error correction rate
             lag.incar.rate+
             #transient terms
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.Females,
           data=ecm.dat)


coeftest(model2.med, vcov = vcovHC(model2.med, type="HC1"))

set.seed(21093)
females.med<-mediate(model.m=model2.med,
                    model.y=model2,
                    treat="lag.punitiveness.Females",
                    mediator="lag.punitiveness.Males",
                    long=FALSE,
                    robustSE=T,sims=10000)











###run model
ecm.dat$lag.punitiveness.65<-ecm.dat$`lag.punitiveness.65+`


model13<-lm(diff.incar.rate~
                 #error correction rate
                 lag.incar.rate+
                 #transient terms
                 diff.Gini+
                 diff.prct.repub+
                 diff.viol.crime.rate+
                 #equlibrium terms
                 lag.pct.black+
                 lag.mood+
                 lag.Gini+
                 lag.prct.repub+
                 lag.unemp.rate+
                 lag.viol.crime.rate+
                 lag.punitiveness.18+
                 lag.punitiveness.45+
            #     lag.punitiveness.30+
                 lag.punitiveness.65,
               data=ecm.dat)

coeftest(model13, vcov = vcovHC(model13, type="HC1"))



model14<-lm(diff.incar.rate~
              #error correction rate
              lag.incar.rate+
              #transient terms
              diff.Gini+
              diff.prct.repub+
              diff.viol.crime.rate+
              #equlibrium terms
              lag.pct.black+
              lag.mood+
              lag.Gini+
              lag.prct.repub+
              lag.unemp.rate+
              lag.viol.crime.rate+
              #lag.punitiveness.18+
              #lag.punitiveness.45+
              lag.punitiveness.30,#+
              #lag.punitiveness.65,
            data=ecm.dat)

coeftest(model14, vcov = vcovHC(model14, type="HC1"))


model3.med<-lm(lag.punitiveness.30~
             #error correction rate
             lag.incar.rate+
             #transient terms
             diff.Gini+
             diff.prct.repub+
             diff.viol.crime.rate+
             #equlibrium terms
             lag.pct.black+
             lag.mood+
             lag.Gini+
             lag.prct.repub+
             lag.unemp.rate+
             lag.viol.crime.rate+
             lag.punitiveness.18+
             lag.punitiveness.45+
             lag.punitiveness.65,
           data=ecm.dat)



coeftest(model3.med, vcov = vcovHC(model3.med, type="HC1"))

set.seed(21093)
age.18.med<-mediate(model.m=model3.med,
                     model.y=model3,
                     treat="lag.punitiveness.18",
                     mediator="lag.punitiveness.30",
                     long=FALSE,
                     robustSE=T,sims=10000)


set.seed(21093)
age.45.med<-mediate(model.m=model3.med,
                    model.y=model3,
                    treat="lag.punitiveness.45",
                    mediator="lag.punitiveness.30",
                    long=FALSE,
                    robustSE=T,sims=10000)




set.seed(21093)
age.65.med<-mediate(model.m=model3.med,
                    model.y=model3,
                    treat="lag.punitiveness.65",
                    mediator="lag.punitiveness.30",
                    long=FALSE,
                    robustSE=T,sims=10000)







#VAR model
###run model
model.var<-lm(diff.incar.rate~
             #error correction rate
             lag.incar.rate,
           data=ecm.dat)






###correlation by group
require(plyr)
cor.by.group <- function(punitiveness.by.group)
{
  return(data.frame(COR = cor(punitiveness.by.group$punitiveness.blacks, punitiveness.by.group$punitiveness.whites)))
}

race.out<-ddply(punitiveness.by.group, .(fipsst), cor.by.group)




require(plyr)
cor.by.group <- function(punitiveness.by.group)
{
  return(data.frame(COR = cor(punitiveness.by.group$punitiveness.Males, 
                              punitiveness.by.group$punitiveness.Females)))
}

sex.out<-ddply(punitiveness.by.group, .(fipsst), cor.by.group)




require(plyr)
cor.by.group <- function(punitiveness.by.group)
{
  return(data.frame(COR = cor(punitiveness.by.group$punitiveness.18, 
                              punitiveness.by.group$punitiveness.30)))
}

age.out1<-ddply(punitiveness.by.group, .(fipsst), cor.by.group)






require(plyr)
cor.by.group <- function(punitiveness.by.group)
{
  return(data.frame(COR = cor(punitiveness.by.group$punitiveness.18, 
                              punitiveness.by.group$punitiveness.45)))
}

age.out2<-ddply(punitiveness.by.group, .(fipsst), cor.by.group)





require(plyr)
cor.by.group <- function(punitiveness.by.group)
{
  return(data.frame(COR = cor(punitiveness.by.group$punitiveness.18, 
                              punitiveness.by.group$`punitiveness.65+`)))
}

age.out3<-ddply(punitiveness.by.group, .(fipsst), cor.by.group)





require(plyr)
cor.by.group <- function(punitiveness.by.group)
{
  return(data.frame(COR = cor(punitiveness.by.group$punitiveness.30, 
                              punitiveness.by.group$punitiveness.45)))
}

age.out4<-ddply(punitiveness.by.group, .(fipsst), cor.by.group)






require(plyr)
cor.by.group <- function(punitiveness.by.group)
{
  return(data.frame(COR = cor(punitiveness.by.group$punitiveness.30, 
                              punitiveness.by.group$`punitiveness.65+`)))
}

age.out5<-ddply(punitiveness.by.group, .(fipsst), cor.by.group)





require(plyr)
cor.by.group <- function(punitiveness.by.group)
{
  return(data.frame(COR = cor(punitiveness.by.group$punitiveness.45, 
                              punitiveness.by.group$`punitiveness.65+`)))
}

age.out6<-ddply(punitiveness.by.group, .(fipsst), cor.by.group)






age.out<-rbind(age.out1,age.out2,age.out3,age.out4,age.out5,age.out6)

